Remove contaminants and filter samples with less than 1000 reads
Paths and libraries setting
# Load main packages, paths and custom functions
source("../../source/main_packages.R")
source("../../source/functions.R")
# Load supplementary packages
packages <- c("decontam", "kableExtra", "microbiome", "cowplot")
invisible(lapply(packages, require, character.only = TRUE))
# create output folders if needed
if (!dir.exists("../../../output/1_MED/1D")) {dir.create("../../../output/1_MED/1D")}Import MED phyloseq object
ps <- readRDS("../../../output/1_MED/1C/1C_MED_phyloseq.rds")Run Decontam
Detect contaminants
# set threshold
threshold=0.6
# Preprocess
df <- as.data.frame(sample_data(ps)) # Put sample_data into a ggplot-friendly data.frame
df <- data.frame(sample_data(ps)) # Put sample_data into a ggplot-friendly data.frame
df <- df[order(df$Read_depth),]
df$Index <- seq(nrow(df))
# read depth
p <- ggplot(data=df, aes(x=Index, y=Read_depth, color=Control)) + geom_point()
p # Identify contaminants - Prevalence method
sample_data(ps)$is.neg <- sample_data(ps)$Control == "Control sample"
contamdf.prev <- isContaminant(ps, method="prevalence", neg="is.neg", threshold=threshold)
# table of positive and false contaminant
table(contamdf.prev$contaminant)##
## FALSE TRUE
## 67 3
# contaminant nodes
decontam_asv_MED <- row.names(contamdf.prev[contamdf.prev$contaminant == TRUE, ])
# make phyloseq object of presence-absence in negative controls and true samples
ps.pa <- transform_sample_counts(ps, function(abund) 1*(abund>0))
ps.pa.neg <- prune_samples(sample_data(ps.pa)$Control == "Control sample", ps.pa)
ps.pa.pos <- prune_samples(sample_data(ps.pa)$Control == "True sample", ps.pa)
# make dataframe from phyloseq objects of presence-absence
df.pa <- data.frame(pa.pos=taxa_sums(ps.pa.pos), pa.neg=taxa_sums(ps.pa.neg),
contaminant=contamdf.prev$contaminant)
# plot
p2 <- ggplot(data=df.pa, aes(x=pa.neg, y=pa.pos, color=contaminant)) + geom_point() +
xlab("Prevalence (Negative Controls)") + ylab("Prevalence (True Samples)")
p2# write the contaminant ASV table detected with prevalence
tax <- tax_table(ps) %>% as.matrix()
decontam_df_MED <- tax[row.names(tax) %in% decontam_asv_MED , ]
# print contaminants
decontam_df_MED %>%
kbl() %>%
kable_paper("hover", full_width = F)| Kingdom | Phylum | Class | Order | Family | Genus | Species | |
|---|---|---|---|---|---|---|---|
| N0005 | Bacteria | Proteobacteria | Gammaproteobacteria | Pseudomonadales | Moraxellaceae | Enhydrobacter | NA |
| N0315 | Bacteria | Proteobacteria | Gammaproteobacteria | Enterobacterales | Yersiniaceae | Rahnella1 | NA |
| N0852 | Bacteria | Proteobacteria | Gammaproteobacteria | Pseudomonadales | Pseudomonadaceae | Pseudomonas | NA |
# histogram of the scores
p <- hist(contamdf.prev$p, 100, ylim = c(0,25), xlim = c(0,1), main="", xlab="p", ylab="Frequency")Save contaminant table
Remove contaminants from the phyloseq object
# remove contaminants ASV
alltaxa <- taxa_names(ps)
decontam_taxa <- alltaxa[!(alltaxa %in% decontam_asv_MED)]
ps2 <- prune_taxa(decontam_taxa, ps)
# check ps
ps2 <- check_ps(ps2)
ps2## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 67 taxa and 123 samples ]
## sample_data() Sample Data: [ 123 samples by 12 sample variables ]
## tax_table() Taxonomy Table: [ 67 taxa by 7 taxonomic ranks ]
## refseq() DNAStringSet: [ 67 reference sequences ]
compute_read_counts(ps2)## [1] 6469784
compute_read_counts(ps)-compute_read_counts(ps2)## [1] 17460
ps <- ps2Remove blanks
# blanks
ps.blanks <- subset_samples(ps, Strain=="Blank")
ps.blanks <- check_ps(ps.blanks)
ps.blanks## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 24 taxa and 8 samples ]
## sample_data() Sample Data: [ 8 samples by 12 sample variables ]
## tax_table() Taxonomy Table: [ 24 taxa by 7 taxonomic ranks ]
## refseq() DNAStringSet: [ 24 reference sequences ]
compute_read_counts(ps.blanks)## [1] 16211
# supprimer blanks
ps.filter <- subset_samples(ps, Strain!="Blank")
ps.filter <- check_ps(ps.filter)
ps.filter## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 67 taxa and 115 samples ]
## sample_data() Sample Data: [ 115 samples by 12 sample variables ]
## tax_table() Taxonomy Table: [ 67 taxa by 7 taxonomic ranks ]
## refseq() DNAStringSet: [ 67 reference sequences ]
compute_read_counts(ps.filter)## [1] 6453573
# check nreads
summarize_phyloseq(ps.filter)## Compositional = NO2
## 1] Min. number of reads = 952] Max. number of reads = 7361593] Total number of reads = 64535734] Average number of reads = 56118.02608695655] Median number of reads = 251717] Sparsity = 0.6486696950032456] Any OTU sum to 1 or less? NO8] Number of singletons = 09] Percent of OTUs that are singletons
## (i.e. exactly one read detected across all samples)010] Number of sample variables are: 12SampleStrainFieldCountryOrganSpeciesRunControlSpecies_italicStrain_italicRead_depthis.neg2
## [[1]]
## [1] "1] Min. number of reads = 95"
##
## [[2]]
## [1] "2] Max. number of reads = 736159"
##
## [[3]]
## [1] "3] Total number of reads = 6453573"
##
## [[4]]
## [1] "4] Average number of reads = 56118.0260869565"
##
## [[5]]
## [1] "5] Median number of reads = 25171"
##
## [[6]]
## [1] "7] Sparsity = 0.648669695003245"
##
## [[7]]
## [1] "6] Any OTU sum to 1 or less? NO"
##
## [[8]]
## [1] "8] Number of singletons = 0"
##
## [[9]]
## [1] "9] Percent of OTUs that are singletons \n (i.e. exactly one read detected across all samples)0"
##
## [[10]]
## [1] "10] Number of sample variables are: 12"
##
## [[11]]
## [1] "Sample" "Strain" "Field" "Country"
## [5] "Organ" "Species" "Run" "Control"
## [9] "Species_italic" "Strain_italic" "Read_depth" "is.neg"
Counts
Create dataframe
x <- c("Culex pipiens", "Field - Bosc", "Field - Camping Europe", "Laboratory - Lavar",
"Culex quinquefasciatus", "Field - Guadeloupe", "Laboratory - Slab TC (Wolbachia -)",
"Aedes aegypti (Guadeloupe)",
"Total")
y <- c("Reads", "Oligotypes", "Samples")
df <- data.frame(matrix(ncol=3, nrow=9))
rownames(df) <- x
colnames(df) <- y
df2 <- df
df3 <- dfWhole + ovary
# Culex pipiens
ps.pipiens <- subset_samples(ps.filter, Species=="Culex pipiens")
ps.pipiens <- check_ps(ps.pipiens)
## All Strain
### oligotype
nrow(ps.pipiens@otu_table) -> df["Culex pipiens", "Oligotypes"]
### samples
nrow(ps.pipiens@sam_data) -> df["Culex pipiens", "Samples"]
### reads
compute_read_counts(ps.pipiens) -> df["Culex pipiens", "Reads"]
## Bosc
### oligotype
nrow(otu_table(ps.pipiens %>%
subset_samples(Strain=="Field - Bosc") %>%
check_ps())) -> df["Field - Bosc", "Oligotypes"]
### samples
nrow(ps.pipiens@sam_data %>% subset_samples(Strain=="Field - Bosc")) -> df["Field - Bosc", "Samples"]
### reads
compute_read_counts(ps.pipiens %>% subset_samples(Strain=="Field - Bosc")) -> df["Field - Bosc", "Reads"]
## Camping Europe
### oligotype
nrow(otu_table(ps.pipiens %>%
subset_samples(Strain=="Field - Camping Europe") %>%
check_ps())) -> df["Field - Camping Europe", "Oligotypes"]
### samples
nrow(ps.pipiens@sam_data %>% subset_samples(Strain=="Field - Camping Europe")) -> df["Field - Camping Europe", "Samples"]
### reads
compute_read_counts(ps.pipiens %>% subset_samples(Strain=="Field - Camping Europe")) -> df["Field - Camping Europe", "Reads"]
## Lavar (labo)
### oligotype
nrow(otu_table(ps.pipiens %>%
subset_samples(Strain=="Laboratory - Lavar") %>%
check_ps())) -> df["Laboratory - Lavar", "Oligotypes"]
### samples
nrow(ps.pipiens@sam_data %>% subset_samples(Strain=="Laboratory - Lavar")) -> df["Laboratory - Lavar", "Samples"]
### reads
compute_read_counts(ps.pipiens %>% subset_samples(Strain=="Laboratory - Lavar")) -> df["Laboratory - Lavar", "Reads"]
# Culex quinquefasciatus
ps.quinque <- subset_samples(ps.filter, Species=="Culex quinquefasciatus")
ps.quinque <- check_ps(ps.quinque)
## All Strain
### oligotype
nrow(ps.quinque@otu_table) -> df["Culex quinquefasciatus", "Oligotypes"]
### samples
nrow(ps.quinque@sam_data) -> df["Culex quinquefasciatus", "Samples"]
### reads
compute_read_counts(ps.quinque) -> df["Culex quinquefasciatus", "Reads"]
## Guadeloupe
### oligotype
nrow(otu_table(ps.quinque %>%
subset_samples(Strain=="Field - Guadeloupe") %>%
check_ps())) -> df["Field - Guadeloupe", "Oligotypes"]
### samples
nrow(ps.quinque@sam_data %>% subset_samples(Strain=="Field - Guadeloupe")) -> df["Field - Guadeloupe", "Samples"]
### reads
compute_read_counts(ps.quinque %>% subset_samples(Strain=="Field - Guadeloupe")) -> df["Field - Guadeloupe", "Reads"]
## Wolbachia -
### oligotype
nrow(otu_table(ps.quinque %>%
subset_samples(Strain=="Laboratory - Slab TC (Wolbachia -)") %>%
check_ps())) -> df["Laboratory - Slab TC (Wolbachia -)", "Oligotypes"]
### samples
nrow(ps.quinque@sam_data %>% subset_samples(Strain=="Laboratory - Slab TC (Wolbachia -)")) -> df["Laboratory - Slab TC (Wolbachia -)", "Samples"]
### reads
compute_read_counts(ps.quinque %>% subset_samples(Strain=="Laboratory - Slab TC (Wolbachia -)")) -> df["Laboratory - Slab TC (Wolbachia -)", "Reads"]
# Aedes aegypti
ps.aedes <- subset_samples(ps.filter, Species=="Aedes aegypti")
ps.aedes <- check_ps(ps.aedes)
## Guadeloupe
### oligotype
nrow(otu_table(ps.aedes %>%
subset_samples(Strain=="Field - Guadeloupe") %>%
check_ps())) -> df["Aedes aegypti (Guadeloupe)", "Oligotypes"]
### samples
nrow(ps.aedes@sam_data %>% subset_samples(Strain=="Field - Guadeloupe")) -> df["Aedes aegypti (Guadeloupe)", "Samples"]
### reads
compute_read_counts(ps.aedes %>% subset_samples(Strain=="Field - Guadeloupe")) -> df["Aedes aegypti (Guadeloupe)", "Reads"]
# Total
### oligotype
nrow(ps.filter@otu_table) -> df["Total", "Oligotypes"]
### samples
nrow(ps.filter@sam_data) -> df["Total", "Samples"]
### reads
compute_read_counts(ps.filter) -> df["Total", "Reads"]
df %>%
kbl() %>%
kable_paper("hover", full_width = F)| Reads | Oligotypes | Samples | |
|---|---|---|---|
| Culex pipiens | 2574603 | 49 | 83 |
| Field - Bosc | 803574 | 36 | 23 |
| Field - Camping Europe | 869617 | 39 | 25 |
| Laboratory - Lavar | 901412 | 48 | 35 |
| Culex quinquefasciatus | 1925054 | 59 | 21 |
| Field - Guadeloupe | 1721053 | 47 | 7 |
| Laboratory - Slab TC (Wolbachia -) | 204001 | 45 | 14 |
| Aedes aegypti (Guadeloupe) | 1953916 | 54 | 11 |
| Total | 6453573 | 67 | 115 |
Whole
# Culex pipiens
ps.pipiens.whole <- ps.pipiens %>% subset_samples(Organ=="Whole")
ps.pipiens.whole <- ps.pipiens.whole %>% check_ps()
## All Strain
### oligotype
nrow(ps.pipiens.whole@otu_table) -> df2["Culex pipiens", "Oligotypes"]
### samples
nrow(ps.pipiens.whole@sam_data) -> df2["Culex pipiens", "Samples"]
### reads
compute_read_counts(ps.pipiens.whole) -> df2["Culex pipiens", "Reads"]
## Bosc
### oligotype
nrow(otu_table(ps.pipiens.whole %>%
subset_samples(Strain=="Field - Bosc") %>%
check_ps())) -> df2["Field - Bosc", "Oligotypes"]
### samples
nrow(ps.pipiens.whole@sam_data %>% subset_samples(Strain=="Field - Bosc")) -> df2["Field - Bosc", "Samples"]
### reads
compute_read_counts(ps.pipiens.whole %>% subset_samples(Strain=="Field - Bosc")) -> df2["Field - Bosc", "Reads"]
## Camping Europe
### oligotype
nrow(otu_table(ps.pipiens.whole %>%
subset_samples(Strain=="Field - Camping Europe") %>%
check_ps())) -> df2["Field - Camping Europe", "Oligotypes"]
### samples
nrow(ps.pipiens.whole@sam_data %>% subset_samples(Strain=="Field - Camping Europe")) -> df2["Field - Camping Europe", "Samples"]
### reads
compute_read_counts(ps.pipiens.whole %>% subset_samples(Strain=="Field - Camping Europe")) -> df2["Field - Camping Europe", "Reads"]
## Lavar (labo)
### oligotype
nrow(otu_table(ps.pipiens.whole %>%
subset_samples(Strain=="Laboratory - Lavar") %>%
check_ps())) -> df2["Laboratory - Lavar", "Oligotypes"]
### samples
nrow(ps.pipiens.whole@sam_data %>% subset_samples(Strain=="Laboratory - Lavar")) -> df2["Laboratory - Lavar", "Samples"]
### reads
compute_read_counts(ps.pipiens.whole %>% subset_samples(Strain=="Laboratory - Lavar")) -> df2["Laboratory - Lavar", "Reads"]
# Culex quinquefasciatus
ps.quinque.whole <- subset_samples(ps.quinque, Organ=="Whole")
ps.quinque.whole <- check_ps(ps.quinque.whole)
## All Strain
### oligotype
nrow(ps.quinque.whole@otu_table) -> df2["Culex quinquefasciatus", "Oligotypes"]
### samples
nrow(ps.quinque.whole@sam_data) -> df2["Culex quinquefasciatus", "Samples"]
### reads
compute_read_counts(ps.quinque.whole) -> df2["Culex quinquefasciatus", "Reads"]
## Guadeloupe
### oligotype
nrow(otu_table(ps.quinque.whole %>%
subset_samples(Strain=="Field - Guadeloupe") %>%
check_ps())) -> df2["Field - Guadeloupe", "Oligotypes"]
### samples
nrow(ps.quinque.whole@sam_data %>% subset_samples(Strain=="Field - Guadeloupe")) -> df2["Field - Guadeloupe", "Samples"]
### reads
compute_read_counts(ps.quinque.whole %>% subset_samples(Strain=="Field - Guadeloupe")) -> df2["Field - Guadeloupe", "Reads"]
## Wolbachia -
### oligotype
nrow(otu_table(ps.quinque.whole %>%
subset_samples(Strain=="Laboratory - Slab TC (Wolbachia -)") %>%
check_ps())) -> df2["Laboratory - Slab TC (Wolbachia -)", "Oligotypes"]
### samples
nrow(ps.quinque.whole@sam_data %>% subset_samples(Strain=="Laboratory - Slab TC (Wolbachia -)")) -> df2["Laboratory - Slab TC (Wolbachia -)", "Samples"]
### reads
compute_read_counts(ps.quinque.whole %>% subset_samples(Strain=="Laboratory - Slab TC (Wolbachia -)")) -> df2["Laboratory - Slab TC (Wolbachia -)", "Reads"]
# Aedes aegypti
ps.aedes.whole <- subset_samples(ps.aedes, Organ=="Whole")
ps.aedes.whole <- check_ps(ps.aedes.whole)
## Guadeloupe
### oligotype
nrow(otu_table(ps.aedes.whole %>%
subset_samples(Strain=="Field - Guadeloupe") %>%
check_ps())) -> df2["Aedes aegypti (Guadeloupe)", "Oligotypes"]
### samples
nrow(ps.aedes.whole@sam_data %>% subset_samples(Strain=="Field - Guadeloupe")) -> df2["Aedes aegypti (Guadeloupe)", "Samples"]
### reads
compute_read_counts(ps.aedes.whole %>% subset_samples(Strain=="Field - Guadeloupe")) -> df2["Aedes aegypti (Guadeloupe)", "Reads"]
# Total
ps.whole <- ps.filter %>% subset_samples(Organ=="Whole")
ps.whole <- ps.whole %>% check_ps()
### oligotype
nrow(ps.whole@otu_table) -> df2["Total", "Oligotypes"]
### samples
nrow(ps.whole@sam_data) -> df2["Total", "Samples"]
### reads
compute_read_counts(ps.whole) -> df2["Total", "Reads"]
df2 %>%
kbl() %>%
kable_paper("hover", full_width = F)| Reads | Oligotypes | Samples | |
|---|---|---|---|
| Culex pipiens | 1251766 | 45 | 42 |
| Field - Bosc | 545790 | 33 | 13 |
| Field - Camping Europe | 98804 | 34 | 7 |
| Laboratory - Lavar | 607172 | 45 | 22 |
| Culex quinquefasciatus | 205761 | 57 | 18 |
| Field - Guadeloupe | 1760 | 38 | 4 |
| Laboratory - Slab TC (Wolbachia -) | 204001 | 45 | 14 |
| Aedes aegypti (Guadeloupe) | 1945807 | 54 | 9 |
| Total | 3403334 | 67 | 69 |
Ovary
# Culex pipiens
ps.pipiens.ovary <- ps.pipiens %>% subset_samples(Organ=="Ovary")
ps.pipiens.ovary <- ps.pipiens.ovary %>% check_ps()
## All Strain
### oligotype
nrow(ps.pipiens.ovary@otu_table) -> df3["Culex pipiens", "Oligotypes"]
### samples
nrow(ps.pipiens.ovary@sam_data) -> df3["Culex pipiens", "Samples"]
### reads
compute_read_counts(ps.pipiens.ovary) -> df3["Culex pipiens", "Reads"]
## Bosc
### oligotype
nrow(otu_table(ps.pipiens.ovary %>%
subset_samples(Strain=="Field - Bosc") %>%
check_ps())) -> df3["Field - Bosc", "Oligotypes"]
### samples
nrow(ps.pipiens.ovary@sam_data %>% subset_samples(Strain=="Field - Bosc")) -> df3["Field - Bosc", "Samples"]
### reads
compute_read_counts(ps.pipiens.ovary %>% subset_samples(Strain=="Field - Bosc")) -> df3["Field - Bosc", "Reads"]
## Camping Europe
### oligotype
nrow(otu_table(ps.pipiens.ovary %>%
subset_samples(Strain=="Field - Camping Europe") %>%
check_ps())) -> df3["Field - Camping Europe", "Oligotypes"]
### samples
nrow(ps.pipiens.ovary@sam_data %>% subset_samples(Strain=="Field - Camping Europe")) -> df3["Field - Camping Europe", "Samples"]
### reads
compute_read_counts(ps.pipiens.ovary %>% subset_samples(Strain=="Field - Camping Europe")) -> df3["Field - Camping Europe", "Reads"]
## Lavar (labo)
### oligotype
nrow(otu_table(ps.pipiens.ovary %>%
subset_samples(Strain=="Laboratory - Lavar") %>%
check_ps())) -> df3["Laboratory - Lavar", "Oligotypes"]
### samples
nrow(ps.pipiens.ovary@sam_data %>% subset_samples(Strain=="Laboratory - Lavar")) -> df3["Laboratory - Lavar", "Samples"]
### reads
compute_read_counts(ps.pipiens.ovary %>% subset_samples(Strain=="Laboratory - Lavar")) -> df3["Laboratory - Lavar", "Reads"]
# Culex quinquefasciatus
ps.quinque.ovary <- subset_samples(ps.quinque, Organ=="Ovary")
ps.quinque.ovary <- check_ps(ps.quinque.ovary)
## All Strain
### oligotype
nrow(ps.quinque.ovary@otu_table) -> df3["Culex quinquefasciatus", "Oligotypes"]
### samples
nrow(ps.quinque.ovary@sam_data) -> df3["Culex quinquefasciatus", "Samples"]
### reads
compute_read_counts(ps.quinque.ovary) -> df3["Culex quinquefasciatus", "Reads"]
## Guadeloupe
### oligotype
nrow(otu_table(ps.quinque.ovary %>%
subset_samples(Strain=="Field - Guadeloupe") %>%
check_ps())) -> df3["Field - Guadeloupe", "Oligotypes"]
### samples
nrow(ps.quinque.ovary@sam_data %>% subset_samples(Strain=="Field - Guadeloupe")) -> df3["Field - Guadeloupe", "Samples"]
### reads
compute_read_counts(ps.quinque.ovary %>% subset_samples(Strain=="Field - Guadeloupe")) -> df3["Field - Guadeloupe", "Reads"]
## Wolbachia -
### oligotype
# nrow(otu_table(ps.quinque.ovary %>%
# subset_samples(Strain=="Laboratory - Slab TC (Wolbachia -)") %>%
# check_ps())) -> df3["Laboratory - Slab TC (Wolbachia -)", "Oligotypes"]
### samples
# nrow(ps.quinque.ovary@sam_data %>% subset_samples(Strain=="Laboratory - Slab TC (Wolbachia -)")) -> df3["Laboratory - Slab TC (Wolbachia -)", "Samples"]
# ### reads
# compute_read_counts(ps.quinque.ovary %>% subset_samples(Strain=="Laboratory - Slab TC (Wolbachia -)")) -> df3["Laboratory - Slab TC (Wolbachia -)", "Reads"]
# Aedes aegypti
ps.aedes.ovary <- subset_samples(ps.aedes, Organ=="Ovary")
ps.aedes.ovary <- check_ps(ps.aedes.ovary)
## Guadeloupe
### oligotype
nrow(otu_table(ps.aedes.ovary %>%
subset_samples(Strain=="Field - Guadeloupe") %>%
check_ps())) -> df3["Aedes aegypti (Guadeloupe)", "Oligotypes"]
### samples
nrow(ps.aedes.ovary@sam_data %>% subset_samples(Strain=="Field - Guadeloupe")) -> df3["Aedes aegypti (Guadeloupe)", "Samples"]
### reads
compute_read_counts(ps.aedes.ovary %>% subset_samples(Strain=="Field - Guadeloupe")) -> df3["Aedes aegypti (Guadeloupe)", "Reads"]
# Total
ps.ovary <- ps.filter %>% subset_samples(Organ=="Ovary")
ps.ovary <- ps.ovary %>% check_ps()
### oligotype
nrow(ps.ovary@otu_table) -> df3["Total", "Oligotypes"]
### samples
nrow(ps.ovary@sam_data) -> df3["Total", "Samples"]
### reads
compute_read_counts(ps.ovary) -> df3["Total", "Reads"]
df3[is.na(df3)] <- 0
df3 %>%
kbl() %>%
kable_paper("hover", full_width = F)| Reads | Oligotypes | Samples | |
|---|---|---|---|
| Culex pipiens | 1322837 | 34 | 41 |
| Field - Bosc | 257784 | 26 | 10 |
| Field - Camping Europe | 770813 | 27 | 18 |
| Laboratory - Lavar | 294240 | 32 | 13 |
| Culex quinquefasciatus | 1719293 | 26 | 3 |
| Field - Guadeloupe | 1719293 | 26 | 3 |
| Laboratory - Slab TC (Wolbachia -) | 0 | 0 | 0 |
| Aedes aegypti (Guadeloupe) | 8109 | 30 | 2 |
| Total | 3050239 | 48 | 46 |
Plot of number of reads by sample
sample_data(ps.filter)$Read_depth <- sample_sums(ps.filter)
metadata_read <- data.frame(ps.filter@sam_data)
p <- ggplot(metadata_read, aes(x = Sample, y = Read_depth))+
geom_bar(position = "dodge", stat = "identity")+
scale_fill_manual(values = col)+
theme_bw() +
theme(axis.text.x = element_text(angle = 90)) +
ggtitle("") +
theme(legend.title = element_text(size = 18),
legend.position="bottom",
legend.text=element_text(size=11),
panel.spacing.y=unit(1, "lines"),
panel.spacing.x=unit(0.8, "lines"),
panel.spacing=unit(0,"lines"),
strip.background=element_rect(color="grey30", fill="grey90"),
panel.border=element_rect(color="grey90"),
axis.ticks.x=element_blank(),
strip.text.x = element_text(size = 12)) +
facet_wrap(~Species_italic+Strain_italic+Organ, scales = "free_x", ncol=4, labeller=label_parsed)+
labs(y="Sequence counts")+
ylim(0, 550000)+
geom_text(aes(label=Read_depth), position=position_dodge(width=0.5), size=2, hjust=-0.25, vjust=0.25, angle=90)+
guides(fill=guide_legend(title="Oligotype", label.theme = element_text(size = 10, face = "italic", colour = "Black", angle = 0)))
pFilter samples
Remove sample with number of reads <1000
a <- prune_samples(sample_sums(ps.filter)<=1000, ps.filter)
test <- data.frame(a@sam_data)
test <- test[test$Sample!="NP20" & test$Sample!="NP29" & test$Sample!="NP30" & test$Sample!="NP34" & test$Sample!="NP36",]
toremove <- test$Sample
ps.filter <- subset_samples(ps.filter, !(Sample %in% toremove))
ps.filter## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 67 taxa and 113 samples ]
## sample_data() Sample Data: [ 113 samples by 12 sample variables ]
## tax_table() Taxonomy Table: [ 67 taxa by 7 taxonomic ranks ]
## refseq() DNAStringSet: [ 67 reference sequences ]
compute_read_counts(ps.filter)## [1] 6452623
Plot of final number of reads by sample
sample_data(ps.filter)$Read_depth <- sample_sums(ps.filter)
metadata_read <- data.frame(ps.filter@sam_data)
metadata_read_whole <- metadata_read[metadata_read$Organ=="Whole",]
metadata_read_ovary <- metadata_read[metadata_read$Organ=="Ovary",]
guide_italics <- guides(fill = guide_legend(label.theme = element_text(size = 10, face = "italic", colour = "Black", angle = 0)))
p <- ggplot(metadata_read_whole, aes(x = Sample, y = Read_depth))+
geom_bar(position = "dodge", stat = "identity")+
scale_fill_manual(values = col)+
theme_bw() +
theme(axis.text.x = element_text(angle = 45, vjust=1, hjust=1, size=15)) +
ggtitle("") +
guide_italics+
theme(legend.title = element_text(size = 18),
legend.position="bottom",
legend.text=element_text(size=16),
panel.spacing.y=unit(1, "lines"),
panel.spacing.x=unit(0.8, "lines"),
panel.spacing=unit(0,"lines"),
strip.background=element_rect(color="grey30", fill="grey90"),
panel.border=element_rect(color="grey90"),
axis.ticks.x=element_blank(),
strip.text.x = element_text(size = 16),
axis.text = element_text(size = 15)) +
facet_wrap(~Species_italic+Strain_italic+Organ, scales = "free_x", ncol=3, labeller=label_parsed)+
labs(y="Sequence counts")+
ylim(0, 1000000)+
geom_text(aes(label=Read_depth), position=position_dodge(width=0.5), size=4, angle=90, hjust=-0.1, vjust=0.25)+
guides(fill=guide_legend(title="Oligotype"))
pp2 <- ggplot(metadata_read_ovary, aes(x = Sample, y = Read_depth))+
geom_bar(position = "dodge", stat = "identity")+
scale_fill_manual(values = col)+
theme_bw() +
theme(axis.text.x = element_text(angle = 45, vjust=1, hjust=1, size=15)) +
ggtitle("") +
guide_italics+
theme(legend.title = element_text(size = 18),
legend.position="bottom",
legend.text=element_text(size=16),
panel.spacing.y=unit(1, "lines"),
panel.spacing.x=unit(0.8, "lines"),
panel.spacing=unit(0,"lines"),
strip.background=element_rect(color="grey30", fill="grey90"),
panel.border=element_rect(color="grey90"),
axis.ticks.x=element_blank(),
strip.text.x = element_text(size = 16),
axis.text = element_text(size = 15)) +
facet_wrap(~Species_italic+Strain_italic+Organ, scales = "free_x", ncol=3, labeller=label_parsed)+
labs(y="Sequence counts")+
ylim(0, 1000000)+
geom_text(aes(label=Read_depth), position=position_dodge(width=0.5), size=4, angle=90, hjust=-0.1, vjust=0.25)+
guides(fill=guide_legend(title="Oligotype"))
p2# panels
p_group <- plot_grid(p+theme(legend.position="none"),
p2+theme(legend.position="none"),
nrow=2,
ncol=1)+
draw_plot_label(c("A", "B"), c(0, 0), c(1, 0.5), size = 15)
legend_plot <- get_legend(p + theme(legend.position="bottom"))
p_counts <- plot_grid(p_group, legend_plot, nrow=2, ncol=1, rel_heights = c(1, .1))
p_countsCounts after removing samples
Create a dataframe
x <- c("Culex pipiens", "Field - Bosc", "Field - Camping Europe", "Laboratory - Lavar",
"Culex quinquefasciatus", "Field - Guadeloupe", "Laboratory - Slab TC (Wolbachia -)",
"Aedes aegypti (Guadeloupe)",
"Total")
y <- c("Reads", "Oligotypes", "Samples")
df <- data.frame(matrix(ncol=3, nrow=9))
rownames(df) <- x
colnames(df) <- y
df2 <- df
df3 <- dfWhole + Ovary
# Culex pipiens
ps.pipiens <- subset_samples(ps.filter, Species=="Culex pipiens")
ps.pipiens <- check_ps(ps.pipiens)
## All Strain
### oligotype
nrow(ps.pipiens@otu_table) -> df["Culex pipiens", "Oligotypes"]
### samples
nrow(ps.pipiens@sam_data) -> df["Culex pipiens", "Samples"]
### reads
compute_read_counts(ps.pipiens) -> df["Culex pipiens", "Reads"]
## Bosc
### oligotype
nrow(otu_table(ps.pipiens %>%
subset_samples(Strain=="Field - Bosc") %>%
check_ps())) -> df["Field - Bosc", "Oligotypes"]
### samples
nrow(ps.pipiens@sam_data %>% subset_samples(Strain=="Field - Bosc")) -> df["Field - Bosc", "Samples"]
### reads
compute_read_counts(ps.pipiens %>% subset_samples(Strain=="Field - Bosc")) -> df["Field - Bosc", "Reads"]
## Camping Europe
### oligotype
nrow(otu_table(ps.pipiens %>%
subset_samples(Strain=="Field - Camping Europe") %>%
check_ps())) -> df["Field - Camping Europe", "Oligotypes"]
### samples
nrow(ps.pipiens@sam_data %>% subset_samples(Strain=="Field - Camping Europe")) -> df["Field - Camping Europe", "Samples"]
### reads
compute_read_counts(ps.pipiens %>% subset_samples(Strain=="Field - Camping Europe")) -> df["Field - Camping Europe", "Reads"]
## Lavar (labo)
### oligotype
nrow(otu_table(ps.pipiens %>%
subset_samples(Strain=="Laboratory - Lavar") %>%
check_ps())) -> df["Laboratory - Lavar", "Oligotypes"]
### samples
nrow(ps.pipiens@sam_data %>% subset_samples(Strain=="Laboratory - Lavar")) -> df["Laboratory - Lavar", "Samples"]
### reads
compute_read_counts(ps.pipiens %>% subset_samples(Strain=="Laboratory - Lavar")) -> df["Laboratory - Lavar", "Reads"]
# Culex quinquefasciatus
ps.quinque <- subset_samples(ps.filter, Species=="Culex quinquefasciatus")
ps.quinque <- check_ps(ps.quinque)
## All Strain
### oligotype
nrow(ps.quinque@otu_table) -> df["Culex quinquefasciatus", "Oligotypes"]
### samples
nrow(ps.quinque@sam_data) -> df["Culex quinquefasciatus", "Samples"]
### reads
compute_read_counts(ps.quinque) -> df["Culex quinquefasciatus", "Reads"]
## Guadeloupe
### oligotype
nrow(otu_table(ps.quinque %>%
subset_samples(Strain=="Field - Guadeloupe") %>%
check_ps())) -> df["Field - Guadeloupe", "Oligotypes"]
### samples
nrow(ps.quinque@sam_data %>% subset_samples(Strain=="Field - Guadeloupe")) -> df["Field - Guadeloupe", "Samples"]
### reads
compute_read_counts(ps.quinque %>% subset_samples(Strain=="Field - Guadeloupe")) -> df["Field - Guadeloupe", "Reads"]
## Wolbachia -
### oligotype
nrow(otu_table(ps.quinque %>%
subset_samples(Strain=="Laboratory - Slab TC (Wolbachia -)") %>%
check_ps())) -> df["Laboratory - Slab TC (Wolbachia -)", "Oligotypes"]
### samples
nrow(ps.quinque@sam_data %>% subset_samples(Strain=="Laboratory - Slab TC (Wolbachia -)")) -> df["Laboratory - Slab TC (Wolbachia -)", "Samples"]
### reads
compute_read_counts(ps.quinque %>% subset_samples(Strain=="Laboratory - Slab TC (Wolbachia -)")) -> df["Laboratory - Slab TC (Wolbachia -)", "Reads"]
# Aedes aegypti
ps.aedes <- subset_samples(ps.filter, Species=="Aedes aegypti")
ps.aedes <- check_ps(ps.aedes)
## Guadeloupe
### oligotype
nrow(otu_table(ps.aedes %>%
subset_samples(Strain=="Field - Guadeloupe") %>%
check_ps())) -> df["Aedes aegypti (Guadeloupe)", "Oligotypes"]
### samples
nrow(ps.aedes@sam_data %>% subset_samples(Strain=="Field - Guadeloupe")) -> df["Aedes aegypti (Guadeloupe)", "Samples"]
### reads
compute_read_counts(ps.aedes %>% subset_samples(Strain=="Field - Guadeloupe")) -> df["Aedes aegypti (Guadeloupe)", "Reads"]
# Total
### oligotype
nrow(ps.filter@otu_table) -> df["Total", "Oligotypes"]
### samples
nrow(ps.filter@sam_data) -> df["Total", "Samples"]
### reads
compute_read_counts(ps.filter) -> df["Total", "Reads"]
df %>%
kbl() %>%
kable_paper("hover", full_width = F)| Reads | Oligotypes | Samples | |
|---|---|---|---|
| Culex pipiens | 2574050 | 49 | 82 |
| Field - Bosc | 803574 | 36 | 23 |
| Field - Camping Europe | 869064 | 38 | 24 |
| Laboratory - Lavar | 901412 | 48 | 35 |
| Culex quinquefasciatus | 1924657 | 59 | 20 |
| Field - Guadeloupe | 1721053 | 47 | 7 |
| Laboratory - Slab TC (Wolbachia -) | 203604 | 45 | 13 |
| Aedes aegypti (Guadeloupe) | 1953916 | 54 | 11 |
| Total | 6452623 | 67 | 113 |
Whole
# Culex pipiens
ps.pipiens.whole <- ps.pipiens %>% subset_samples(Organ=="Whole")
ps.pipiens.whole <- ps.pipiens.whole %>% check_ps()
## All Strain
### oligotype
nrow(ps.pipiens.whole@otu_table) -> df2["Culex pipiens", "Oligotypes"]
### samples
nrow(ps.pipiens.whole@sam_data) -> df2["Culex pipiens", "Samples"]
### reads
compute_read_counts(ps.pipiens.whole) -> df2["Culex pipiens", "Reads"]
## Bosc
### oligotype
nrow(otu_table(ps.pipiens.whole %>%
subset_samples(Strain=="Field - Bosc") %>%
check_ps())) -> df2["Field - Bosc", "Oligotypes"]
### samples
nrow(ps.pipiens.whole@sam_data %>% subset_samples(Strain=="Field - Bosc")) -> df2["Field - Bosc", "Samples"]
### reads
compute_read_counts(ps.pipiens.whole %>% subset_samples(Strain=="Field - Bosc")) -> df2["Field - Bosc", "Reads"]
## Camping Europe
### oligotype
nrow(otu_table(ps.pipiens.whole %>%
subset_samples(Strain=="Field - Camping Europe") %>%
check_ps())) -> df2["Field - Camping Europe", "Oligotypes"]
### samples
nrow(ps.pipiens.whole@sam_data %>% subset_samples(Strain=="Field - Camping Europe")) -> df2["Field - Camping Europe", "Samples"]
### reads
compute_read_counts(ps.pipiens.whole %>% subset_samples(Strain=="Field - Camping Europe")) -> df2["Field - Camping Europe", "Reads"]
## Lavar (labo)
### oligotype
nrow(otu_table(ps.pipiens.whole %>%
subset_samples(Strain=="Laboratory - Lavar") %>%
check_ps())) -> df2["Laboratory - Lavar", "Oligotypes"]
### samples
nrow(ps.pipiens.whole@sam_data %>% subset_samples(Strain=="Laboratory - Lavar")) -> df2["Laboratory - Lavar", "Samples"]
### reads
compute_read_counts(ps.pipiens.whole %>% subset_samples(Strain=="Laboratory - Lavar")) -> df2["Laboratory - Lavar", "Reads"]
# Culex quinquefasciatus
ps.quinque.whole <- subset_samples(ps.quinque, Organ=="Whole")
ps.quinque.whole <- check_ps(ps.quinque.whole)
## All Strain
### oligotype
nrow(ps.quinque.whole@otu_table) -> df2["Culex quinquefasciatus", "Oligotypes"]
### samples
nrow(ps.quinque.whole@sam_data) -> df2["Culex quinquefasciatus", "Samples"]
### reads
compute_read_counts(ps.quinque.whole) -> df2["Culex quinquefasciatus", "Reads"]
## Guadeloupe
### oligotype
nrow(otu_table(ps.quinque.whole %>%
subset_samples(Strain=="Field - Guadeloupe") %>%
check_ps())) -> df2["Field - Guadeloupe", "Oligotypes"]
### samples
nrow(ps.quinque.whole@sam_data %>% subset_samples(Strain=="Field - Guadeloupe")) -> df2["Field - Guadeloupe", "Samples"]
### reads
compute_read_counts(ps.quinque.whole %>% subset_samples(Strain=="Field - Guadeloupe")) -> df2["Field - Guadeloupe", "Reads"]
## Wolbachia -
### oligotype
nrow(otu_table(ps.quinque.whole %>%
subset_samples(Strain=="Laboratory - Slab TC (Wolbachia -)") %>%
check_ps())) -> df2["Laboratory - Slab TC (Wolbachia -)", "Oligotypes"]
### samples
nrow(ps.quinque.whole@sam_data %>% subset_samples(Strain=="Laboratory - Slab TC (Wolbachia -)")) -> df2["Laboratory - Slab TC (Wolbachia -)", "Samples"]
### reads
compute_read_counts(ps.quinque.whole %>% subset_samples(Strain=="Laboratory - Slab TC (Wolbachia -)")) -> df2["Laboratory - Slab TC (Wolbachia -)", "Reads"]
# Aedes aegypti
ps.aedes.whole <- subset_samples(ps.aedes, Organ=="Whole")
ps.aedes.whole <- check_ps(ps.aedes.whole)
## Guadeloupe
### oligotype
nrow(otu_table(ps.aedes.whole %>%
subset_samples(Strain=="Field - Guadeloupe") %>%
check_ps())) -> df2["Aedes aegypti (Guadeloupe)", "Oligotypes"]
### samples
nrow(ps.aedes.whole@sam_data %>% subset_samples(Strain=="Field - Guadeloupe")) -> df2["Aedes aegypti (Guadeloupe)", "Samples"]
### reads
compute_read_counts(ps.aedes.whole %>% subset_samples(Strain=="Field - Guadeloupe")) -> df2["Aedes aegypti (Guadeloupe)", "Reads"]
# Total
ps.whole <- ps.filter %>% subset_samples(Organ=="Whole")
ps.whole <- ps.whole %>% check_ps()
### oligotype
nrow(ps.whole@otu_table) -> df2["Total", "Oligotypes"]
### samples
nrow(ps.whole@sam_data) -> df2["Total", "Samples"]
### reads
compute_read_counts(ps.whole) -> df2["Total", "Reads"]
df2 %>%
kbl() %>%
kable_paper("hover", full_width = F)| Reads | Oligotypes | Samples | |
|---|---|---|---|
| Culex pipiens | 1251213 | 45 | 41 |
| Field - Bosc | 545790 | 33 | 13 |
| Field - Camping Europe | 98251 | 33 | 6 |
| Laboratory - Lavar | 607172 | 45 | 22 |
| Culex quinquefasciatus | 205364 | 57 | 17 |
| Field - Guadeloupe | 1760 | 38 | 4 |
| Laboratory - Slab TC (Wolbachia -) | 203604 | 45 | 13 |
| Aedes aegypti (Guadeloupe) | 1945807 | 54 | 9 |
| Total | 3402384 | 67 | 67 |
Ovary
# Culex pipiens
ps.pipiens.ovary <- ps.pipiens %>% subset_samples(Organ=="Ovary")
ps.pipiens.ovary <- ps.pipiens.ovary %>% check_ps()
## All Strain
### oligotype
nrow(ps.pipiens.ovary@otu_table) -> df3["Culex pipiens", "Oligotypes"]
### samples
nrow(ps.pipiens.ovary@sam_data) -> df3["Culex pipiens", "Samples"]
### reads
compute_read_counts(ps.pipiens.ovary) -> df3["Culex pipiens", "Reads"]
## Bosc
### oligotype
nrow(otu_table(ps.pipiens.ovary %>%
subset_samples(Strain=="Field - Bosc") %>%
check_ps())) -> df3["Field - Bosc", "Oligotypes"]
### samples
nrow(ps.pipiens.ovary@sam_data %>% subset_samples(Strain=="Field - Bosc")) -> df3["Field - Bosc", "Samples"]
### reads
compute_read_counts(ps.pipiens.ovary %>% subset_samples(Strain=="Field - Bosc")) -> df3["Field - Bosc", "Reads"]
## Camping Europe
### oligotype
nrow(otu_table(ps.pipiens.ovary %>%
subset_samples(Strain=="Field - Camping Europe") %>%
check_ps())) -> df3["Field - Camping Europe", "Oligotypes"]
### samples
nrow(ps.pipiens.ovary@sam_data %>% subset_samples(Strain=="Field - Camping Europe")) -> df3["Field - Camping Europe", "Samples"]
### reads
compute_read_counts(ps.pipiens.ovary %>% subset_samples(Strain=="Field - Camping Europe")) -> df3["Field - Camping Europe", "Reads"]
## Lavar (labo)
### oligotype
nrow(otu_table(ps.pipiens.ovary %>%
subset_samples(Strain=="Laboratory - Lavar") %>%
check_ps())) -> df3["Laboratory - Lavar", "Oligotypes"]
### samples
nrow(ps.pipiens.ovary@sam_data %>% subset_samples(Strain=="Laboratory - Lavar")) -> df3["Laboratory - Lavar", "Samples"]
### reads
compute_read_counts(ps.pipiens.ovary %>% subset_samples(Strain=="Laboratory - Lavar")) -> df3["Laboratory - Lavar", "Reads"]
# Culex quinquefasciatus
ps.quinque.ovary <- subset_samples(ps.quinque, Organ=="Ovary")
ps.quinque.ovary <- check_ps(ps.quinque.ovary)
## All Strain
### oligotype
nrow(ps.quinque.ovary@otu_table) -> df3["Culex quinquefasciatus", "Oligotypes"]
### samples
nrow(ps.quinque.ovary@sam_data) -> df3["Culex quinquefasciatus", "Samples"]
### reads
compute_read_counts(ps.quinque.ovary) -> df3["Culex quinquefasciatus", "Reads"]
## Guadeloupe
### oligotype
nrow(otu_table(ps.quinque.ovary %>%
subset_samples(Strain=="Field - Guadeloupe") %>%
check_ps())) -> df3["Field - Guadeloupe", "Oligotypes"]
### samples
nrow(ps.quinque.ovary@sam_data %>% subset_samples(Strain=="Field - Guadeloupe")) -> df3["Field - Guadeloupe", "Samples"]
### reads
compute_read_counts(ps.quinque.ovary %>% subset_samples(Strain=="Field - Guadeloupe")) -> df3["Field - Guadeloupe", "Reads"]
## Wolbachia -
### oligotype
# nrow(otu_table(ps.quinque.ovary %>%
# subset_samples(Strain=="Laboratory - Slab TC (Wolbachia -)") %>%
# check_ps())) -> df3["Laboratory - Slab TC (Wolbachia -)", "Oligotypes"]
### samples
# nrow(ps.quinque.ovary@sam_data %>% subset_samples(Strain=="Laboratory - Slab TC (Wolbachia -)")) -> df3["Laboratory - Slab TC (Wolbachia -)", "Samples"]
# ### reads
# compute_read_counts(ps.quinque.ovary %>% subset_samples(Strain=="Laboratory - Slab TC (Wolbachia -)")) -> df3["Laboratory - Slab TC (Wolbachia -)", "Reads"]
# Aedes aegypti
ps.aedes.ovary <- subset_samples(ps.aedes, Organ=="Ovary")
ps.aedes.ovary <- check_ps(ps.aedes.ovary)
## Guadeloupe
### oligotype
nrow(otu_table(ps.aedes.ovary %>%
subset_samples(Strain=="Field - Guadeloupe") %>%
check_ps())) -> df3["Aedes aegypti (Guadeloupe)", "Oligotypes"]
### samples
nrow(ps.aedes.ovary@sam_data %>% subset_samples(Strain=="Field - Guadeloupe")) -> df3["Aedes aegypti (Guadeloupe)", "Samples"]
### reads
compute_read_counts(ps.aedes.ovary %>% subset_samples(Strain=="Field - Guadeloupe")) -> df3["Aedes aegypti (Guadeloupe)", "Reads"]
# Total
ps.ovary <- ps.filter %>% subset_samples(Organ=="Ovary")
ps.ovary <- ps.ovary %>% check_ps()
### oligotype
nrow(ps.ovary@otu_table) -> df3["Total", "Oligotypes"]
### samples
nrow(ps.ovary@sam_data) -> df3["Total", "Samples"]
### reads
compute_read_counts(ps.ovary) -> df3["Total", "Reads"]
df3[is.na(df3)] <- 0
df3 %>%
kbl() %>%
kable_paper("hover", full_width = F)| Reads | Oligotypes | Samples | |
|---|---|---|---|
| Culex pipiens | 1322837 | 34 | 41 |
| Field - Bosc | 257784 | 26 | 10 |
| Field - Camping Europe | 770813 | 27 | 18 |
| Laboratory - Lavar | 294240 | 32 | 13 |
| Culex quinquefasciatus | 1719293 | 26 | 3 |
| Field - Guadeloupe | 1719293 | 26 | 3 |
| Laboratory - Slab TC (Wolbachia -) | 0 | 0 | 0 |
| Aedes aegypti (Guadeloupe) | 8109 | 30 | 2 |
| Total | 3050239 | 48 | 46 |
Save data
Save count tables and count plots
write.table(df, "../../../output/1_MED/1D/1D_Counts_all.tsv", sep="\t", row.names = TRUE, col.names=NA)
write.table(df2, "../../../output/1_MED/1D/1D_Counts_whole.tsv", sep="\t", row.names = TRUE, col.names=NA)
write.table(df3, "../../../output/1_MED/1D/1D_Counts_ovary.tsv", sep="\t", row.names = TRUE, col.names=NA)
tiff("../../../output/1_MED/1D/1D_counts_by_sample.tiff", units="in", width=20, height=18, res=300)
p_counts
dev.off()## quartz_off_screen
## 2
png("../../../output/1_MED/1D/1D_counts_by_sample.png", units="in", width=20, height=18, res=300)
p_counts
dev.off()## quartz_off_screen
## 2
Save metadata (for Supplementary Table 1)
count_after_decontam <- data.frame(ps.filter@otu_table)
tax_after_decontam <- data.frame(ps.filter@tax_table)
metadata_after_decontam <- data.frame(ps.filter@sam_data)
compute_read_counts(ps.filter)## [1] 6452623
write.table(count_after_decontam, "../../../output/1_MED/1D/1D_count_supplementary_table1.tsv", sep="\t", row.names = TRUE, col.names=NA)
write.table(tax_after_decontam, "../../../output/1_MED/1D/1D_tax_supplementary_table1.tsv", sep="\t", row.names = TRUE, col.names=NA)Save phyloseq object after decontam
saveRDS(ps.filter, "../../../output/1_MED/1D/1D_MED_phyloseq_decontam.rds")Session info
sessionInfo()## R version 3.6.3 (2020-02-29)
## Platform: x86_64-apple-darwin15.6.0 (64-bit)
## Running under: macOS 10.16
##
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/3.6/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/3.6/Resources/lib/libRlapack.dylib
##
## locale:
## [1] fr_FR.UTF-8/fr_FR.UTF-8/fr_FR.UTF-8/C/fr_FR.UTF-8/fr_FR.UTF-8
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] cowplot_1.1.0 microbiome_1.8.0 kableExtra_1.3.1 decontam_1.6.0
## [5] forcats_0.5.0 stringr_1.4.0 dplyr_1.0.2 purrr_0.3.4
## [9] readr_1.4.0 tidyr_1.1.2 tibble_3.0.4 tidyverse_1.3.0
## [13] ggplot2_3.3.2 phyloseq_1.30.0
##
## loaded via a namespace (and not attached):
## [1] nlme_3.1-149 fs_1.5.0 lubridate_1.7.9
## [4] webshot_0.5.2 httr_1.4.2 tools_3.6.3
## [7] backports_1.1.10 bslib_0.2.4 R6_2.4.1
## [10] vegan_2.5-6 DBI_1.1.0 BiocGenerics_0.32.0
## [13] mgcv_1.8-33 colorspace_2.0-0 permute_0.9-5
## [16] ade4_1.7-15 withr_2.3.0 tidyselect_1.1.0
## [19] compiler_3.6.3 cli_2.1.0 rvest_0.3.6
## [22] Biobase_2.46.0 xml2_1.3.2 labeling_0.4.2
## [25] bookdown_0.22 sass_0.3.1 scales_1.1.1
## [28] digest_0.6.26 rmarkdown_2.7 XVector_0.26.0
## [31] pkgconfig_2.0.3 htmltools_0.5.1.1 highr_0.8
## [34] dbplyr_1.4.4 rlang_0.4.10 readxl_1.3.1
## [37] rstudioapi_0.11 farver_2.0.3 jquerylib_0.1.3
## [40] generics_0.0.2 jsonlite_1.7.1 magrittr_1.5
## [43] biomformat_1.14.0 Matrix_1.2-18 fansi_0.4.1
## [46] Rcpp_1.0.5 munsell_0.5.0 S4Vectors_0.24.4
## [49] Rhdf5lib_1.8.0 ape_5.4-1 lifecycle_0.2.0
## [52] stringi_1.5.3 yaml_2.2.1 MASS_7.3-53
## [55] zlibbioc_1.32.0 Rtsne_0.15 rhdf5_2.30.1
## [58] plyr_1.8.6 grid_3.6.3 blob_1.2.1
## [61] parallel_3.6.3 crayon_1.3.4 lattice_0.20-41
## [64] Biostrings_2.54.0 haven_2.3.1 splines_3.6.3
## [67] multtest_2.42.0 hms_0.5.3 ps_1.4.0
## [70] knitr_1.30 pillar_1.4.6 igraph_1.2.6
## [73] reshape2_1.4.4 codetools_0.2-16 stats4_3.6.3
## [76] reprex_0.3.0 glue_1.4.2 evaluate_0.14
## [79] data.table_1.13.2 modelr_0.1.8 vctrs_0.3.4
## [82] rmdformats_1.0.2 foreach_1.5.1 cellranger_1.1.0
## [85] gtable_0.3.0 assertthat_0.2.1 xfun_0.22
## [88] broom_0.7.2 viridisLite_0.3.0 survival_3.2-7
## [91] iterators_1.0.13 IRanges_2.20.2 cluster_2.1.0
## [94] ellipsis_0.3.1